# This file will be called by "MasterFile.R", which sets "parent_path" for file output.
# This files contains the main codes that output figures for different kappa_x, the spending during crisis that cannot be recovered. 
# ------------------ Plot Specifications ---------------------
margins = c(3.5, 3.5, 1, 1)
plot_slippery_slope_two_economy = TRUE
plot_width = 4.2
plot_height = 3.8

for(iter_case in 1:4){
  # # ------------------ Use baseline calibration parameters. -------------
  load('baseline_calibration.Rdata')
  # NO need to recalibrate. Directly use the default values. 
  source("Model scripts-with recovery kappax.R")  # This script solves the model and provide simulation functions 
  # Setting the key parameter and associated file path
  if(iter_case==1){
    par$kappax = 0.8  #  KEY extra moment: kappax is the loss rate of crisis-time spending. 1-kappax is left over.
    subfigure = "a";
  }else if(iter_case==2){
    par$kappax = 0.4
    subfigure = "b";
  }else if(iter_case==3){
    par$kappax = 0.2
    subfigure = "c";
  }else{
    par$kappax = 0.1
    subfigure = "d";
  }
  print(paste0("------------ Master Iteration : kappax = ", par$kappax, " ---------------"))
  # redo the calibration  
  calibration_results=calibration_moments(par)
  equi_sol = calibration_results$equi_sol
  moments = cbind(calibration_results$key_moments, calibration_results$other_moments)
  print("moments:")
  print(moments)
  
  # Then assign values for next processing
  g_bar = par$g_bar; d_bar = par$d_bar; dt = equi_sol$dt;
  AH = par$AH;  AL = par$AL
  lambda=par$lambda; beta=par$beta; r = par$r; eta_H=par$eta_H; eta_L=par$eta_L
  w_grid = equi_sol$w_grid
  f_list = calibration_results$f_list
  solve_individual_optimal_faster = f_list$solve_individual_optimal_faster; F_fun = f_list$F_fun; Delta_output_fun = f_list$Delta_output_fun; solve_equilibrium=f_list$solve_equilibrium; Simulate_model=f_list$Simulate_model; solve_other_moments=f_list$solve_other_moments; 
  muK_fun=f_list$muK_fun; DeltaK_fun=f_list$DeltaK_fun; muw_fun=f_list$muw_fun; Deltaw_fun=f_list$Deltaw_fun; investment_cost_fun=f_list$investment_cost_fun;  post_crisis_fraction_of_capital=f_list$post_crisis_fraction_of_capital; output_drop_and_gbar=f_list$output_drop_and_gbar;  
  solve_welfare = f_list$solve_welfare
  simulation_fun=f_list$simulation_fun
  avg_zeta = par$avg_zeta; 
  qL = equi_sol$qL; qH = equi_sol$qH;  
  
  equi_sol_baseline = equi_sol    # This is the baseline equilibrium solution 
  
  # Basic moments
  equi_sol_no_govt = solve_equilibrium(d_bar, g_bar=0)
  output_drop_and_gbar( par$g_bar,  equi_sol$w_avg, qH = equi_sol$qH, qL=equi_sol$qL  )
  output_drop_and_gbar( 0,  equi_sol$w_avg, qH=equi_sol_no_govt$qH, qL=equi_sol_no_govt$qL  )
  Deltaw_fun( equi_sol$kappaH, equi_sol$kappaL, equi_sol$w_avg )
  Deltaw_fun( equi_sol_no_govt$kappaH, equi_sol_no_govt$kappaL, equi_sol$w_avg )
  equi_sol = solve_other_moments(equi_sol_baseline)
  equi_sol$private_financing / equi_sol$avg_productivity
  
  # -------------------- Figure A15: Slippery Slope in the Same Economy --------------------
  # Government intervention and the avg w
  g_bar_vec = c( seq(0, g_bar ,length.out=20), seq(g_bar+0.005, 0.25, length.out=20) )
  index_baseline = which.min( abs(g_bar_vec-0.14)  )
  GDP_drop_target = -0.1
  # intervention and next intervention scale. 
  govt_funding_next_crisis_vec = rep(0, length(g_bar_vec))
  equi_sol_no_govt = solve_equilibrium(d_bar, g_bar=0)
  equi_sol_no_govt = solve_other_moments(equi_sol_no_govt)
  w0 = equi_sol_no_govt$w_avg
  K0 = 1
  # Think about the actual amount of intervention. 
  govt_credit_provision_vec = g_bar_vec
  govt_credit_provision_next_crisis = g_bar_vec
  w_next_crisis_vec = g_bar_vec
  K_next_crisis_vec = g_bar_vec
  productivity_vec = rep(0,length(g_bar_vec))
  productivity_next_crisis_vec = rep(0,length(g_bar_vec))
  for(iter_g_bar in 1:length(g_bar_vec)){
    print( paste0("Progress:", iter_g_bar/length(g_bar_vec)*100, "%") )
    resultH = post_crisis_fraction_of_capital(equi_sol_no_govt$qH, d_bar, g_bar_vec[iter_g_bar], type="H")
    resultL = post_crisis_fraction_of_capital(equi_sol_no_govt$qL, d_bar, g_bar_vec[iter_g_bar], type="L")
    kappaH = resultH$kappa
    kappaL = resultL$kappa
    govt_credit_provision_vec[iter_g_bar] = w0*resultH$govt_sector_financing + (1-w0)*resultL$govt_sector_financing
    productivity_vec[iter_g_bar] = w0*AH + (1-w0)*AL
    w = w0  + Deltaw_fun(kappaH, kappaL, w0) # immediately after the crisis
    K = K0 * (1 + DeltaK_fun(kappaH, kappaL, w0) )
    for(iter in 1:10/dt){
      # 10 years
      w = w + muw_fun(equi_sol_no_govt$iH, equi_sol_no_govt$iL, w)*dt
      K = K * (1 + muK_fun(equi_sol_no_govt$iH, equi_sol_no_govt$iL, w)*dt )
    }  
    w_next_crisis_vec[iter_g_bar] = w
    K_next_crisis_vec[iter_g_bar] = K
    result = uniroot( function(g_bar) output_drop_and_gbar(g_bar, w, equi_sol_no_govt$qH, equi_sol_no_govt$qL) - GDP_drop_target, interval = c(0, 2)  )
    govt_funding_next_crisis_vec[iter_g_bar] = result$root
    
    resultH= post_crisis_fraction_of_capital(equi_sol_no_govt$qH, d_bar, govt_funding_next_crisis_vec[iter_g_bar], type="H")
    resultL= post_crisis_fraction_of_capital(equi_sol_no_govt$qL, d_bar, govt_funding_next_crisis_vec[iter_g_bar], type="L")
    # Total credit provision
    govt_credit_provision_next_crisis[iter_g_bar] = w*resultH$govt_sector_financing + (1-w)*resultL$govt_sector_financing
    productivity_next_crisis_vec[iter_g_bar] = w*AH + (1-w)*AL
  }
  sensitivity = ( govt_funding_next_crisis_vec[index_baseline] - govt_funding_next_crisis_vec[1] ) / ( g_bar_vec[index_baseline] - g_bar_vec[1] )
  pdf(paste0(FigurePath,"/FigureA15", subfigure, "_same_economy_slippery_slope_credit.pdf"), width=plot_width,height=plot_height, family = "serif")
  par(mfrow=c(1,1), mar=margins)
  MyPlot( g_bar_vec, govt_funding_next_crisis_vec, xlab=TeX("current govt program size $\\bar{g}$"), ylab=TeX("next crisis (in 10 years) govt intervention $\\bar{g}'$") )
  dev.off()
  
  # ----- Figure A16: Slippery Slope Across Two Economies ----------
  w_starting_vec = equi_sol_no_govt$w_steady 
  years_sim = 20
  GDP_drop_target = -0.1
  # Set expectation as non-government intervention!!! 
  for(w_starting in w_starting_vec){
    dt = 1/4
    # dN_vec = c( rep(0,10/dt), 1, rep(0,(years_sim-10)/dt-1) )
    dN_vec = c( 1, rep(0,10/dt), 1, rep(0,(years_sim-10)/dt-2) )
    # dN_vec = rep( 0, years_sim/dt )
    N = first(which(dN_vec==1))
    equi_sol_no_govt = solve_equilibrium( d_bar, g_bar=0 )
    # par1: the parameters for case 1, i.e., expectation is with intervention and there is actually government intervention. The intervention gamma is given by low gamma (lenient intervention)
    equi_sol_with_govt = solve_equilibrium( d_bar, g_bar )
    equi_sol_with_govt$qH = equi_sol_no_govt$qH
    equi_sol_with_govt$qL = equi_sol_no_govt$qL
    equi_sol_with_govt$iH = equi_sol_no_govt$iH
    equi_sol_with_govt$iL = equi_sol_no_govt$iL
    
    simulation_case1 = simulation_fun( equi_sol_with_govt, dN_vec = dN_vec  , dt=dt , years_sim = years_sim, w_init=w_starting)
    simulation_case1[  , actual_govt_intervention := 0 ]  # scale of govt intervention
    simulation_case1[ which(dN_vec==1) , actual_govt_intervention := equi_sol_with_govt$total_borrowing_H*w + equi_sol_with_govt$total_borrowing_L*(1-w)  ]  # scale of govt intervention
    simulation_case1[  , gbar_needed := 0  ]
    simulation_case1[  , total_govt_financing := 0  ]
    simulation_case1[  , total_private_financing := 0  ]
    # par2: the parameters for case 2, i.e., expectation is laissez-faire, and there is laissez-faire at first crisis, but with intervention at the second crisis to achieve the same GDP decline
    simulation_case2 = simulation_fun( equi_sol_no_govt, dN_vec = dN_vec, dt=dt , years_sim = years_sim, w_init=w_starting)
    # Then solve for the gamma that achieves the same GDP decline.
    simulation_case2[  , gbar_needed := 0  ]
    simulation_case2[  , total_govt_financing := 0  ]
    simulation_case2[  , total_private_financing := 0  ]
    
    MyPlot( zeta_vec, equi_sol_with_govt$result_qL$l_fun(zeta_vec) )
    
    for(iter in 1:nrow(simulation_case2)){
      print( paste0( "Progress:", round( iter/nrow(simulation_case2)*100, 0 )  , "%" ) )
      w = simulation_case1[iter]$w
      result = uniroot( function(g_bar) output_drop_and_gbar(g_bar, w, equi_sol_no_govt$qH, equi_sol_no_govt$qL) - GDP_drop_target, interval = c(0, 2)  )
      g_bar_needed = result$root
      simulation_case1[iter]$gbar_needed = g_bar_needed
      resultH = post_crisis_fraction_of_capital( equi_sol_no_govt$qH, equi_sol_no_govt$d_bar, g_bar_needed, type="H"  )  # Then find the total governmetn credit provision
      resultL = post_crisis_fraction_of_capital( equi_sol_no_govt$qL, equi_sol_no_govt$d_bar, g_bar_needed, type="L"  )
      simulation_case1[iter]$total_govt_financing = w*resultH$govt_sector_financing + (1-w)*resultL$govt_sector_financing
      simulation_case1[iter]$total_private_financing = w*resultH$private_sector_financing + (1-w)*resultL$private_sector_financing
      
      w = simulation_case2[iter]$w
      result = uniroot( function(g_bar) output_drop_and_gbar(g_bar, w, equi_sol_no_govt$qH, equi_sol_no_govt$qL) - GDP_drop_target, interval = c(0, 2), tol=10^(-12)  )
      g_bar_needed = result$root
      simulation_case2[iter]$gbar_needed = g_bar_needed
      resultH = post_crisis_fraction_of_capital( equi_sol_no_govt$qH, equi_sol_no_govt$d_bar, g_bar_needed, type="H"  )  # Then find the total governmetn credit provision
      resultL = post_crisis_fraction_of_capital( equi_sol_no_govt$qL, equi_sol_no_govt$d_bar, g_bar_needed, type="L"  )
      simulation_case2[iter]$total_govt_financing = w*resultH$govt_sector_financing + (1-w)*resultL$govt_sector_financing
      simulation_case2[iter]$total_private_financing = w*resultH$private_sector_financing + (1-w)*resultL$private_sector_financing
    }
    pdf(paste0(FigurePath,"/FigureA16", subfigure ,"_two_economy_slippery_slope.pdf"), width=plot_width, height=plot_height, family="serif")
    par(mfrow=c(1,1), mar=margins)
    MyPlot( simulation_case1$T, 100*(simulation_case1$gbar_needed-simulation_case2$gbar_needed )/simulation_case2$gbar_needed, xlab="simulation years", ylab="extra govt intervention needed (%)" )
    dev.off()
  }

  
}


